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Abstract 

Fluctuations of the instantaneous local Lagrangian strain ejj(r,t), measured 



(N 
> 

m 

\Q ' with respect to a static "reference" lattice, are used to obtain accurate esti 

o' 

^^ ' mates of the elastic constants of model solids from atomistic computer simu 

^.^ ' lations. The measured strains are systematically coarse - grained by averag- 

d • ing them within subsystems (of size -Lfe) of a system (of total size L) in the 

canonical ensemble. Using a simple finite size scaling theory we predict the 

is ' behaviour of the fluctuations < eije^i > as a function of Li,/L and extract 

O ■ 

elastic constants of the system in the thermodynamic limit at nonzero tem- 

^ ' perature. Our method is simple to implement, efficient and general enough to 

be able to handle a wide class of model systems including those with singular 
potentials without any essential modification. We illustrate the technique by 
computing isothermal elastic constants of "hard" and "soft" disk triangular 
solids in two dimensions from Monte Carlo and molecular dynamics simula- 
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tions. We compare our results with those from ear her simulations and theory. 
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I. INTRODUCTION 

One is often interested in long length scale and long time scale phenomena in solids ( eg. 
late stage kinetics of solid state phase transformations [|l],0 ; motion of domain walls inter- 
faces 0; fracture 0; friction |^ etc.). Such phenomena are usually described by continuum 
theories. Microscopic simulations ||^ of finite systems, on the other hand, like molecular dy- 
namics, lattice Boltzmann or Monte Carlo, deal with microscopic variables like the positions 
and velocities of constituent particles and together with detailed knowledge of interatomic 
potentials, hope to build up a description of the macro system from a knowledge of these 
micro variables. How does one recover continuum physics from simulating the dynamics of 
A^ particles? This requires a "coarse-graining" procedure in space (for equilibrium) or both 
space and time for non-equilibrium continuum theories. Over what coarse graining length 
and time scale does one recover results consistent with continuum theories? In this paper, we 
attempt to answer these questions for the simplest nontrivial case, namely, a crystalline solid, 
(without any point, line or surface defects [§]) in equilibrium, at a non zero temperature 
far away from phase transitions. We show that coarse graining of microscopic local strain 
fluctuations obtained from configurations generated in a computer simulation, enables us to 
calculate elastic susceptibilities (compliances) as a function of the coarse graining length. 
Detailed finite size scaling analysis of this data yields finally elastic constants of the solid - 
the essential inputs to continuum elasticity theory 0. The strain field (together with defect 
densities) constitute the coarse grained description of a solid essential for understanding 
solid state phase transitions like structural transitions |[l^J^ and melting [|Tl| -[l3[ . 



Calculation of elastic constants from simulations fall into two categories, viz. they are 
obtained either from, thermal averages of fluctuations of the stress or the strain - the so 



called "fluctuation" methods [14], or from the stress - strain curve as computed from a series 



of simulations [0. Fluctuation methods, though requiring longer runs for accumulating 
statistically significant data, are often preferred because the entire matrix of elastic constants 
can be evaluated in a single run, whereas in the latter method for every elastic constant an 



appropriate strain (or stress) have to be applied. Also mapping out the stress - strain curve 
can be treacherous especially for soft systems where the possibility of setting up a plastic flow 
in the system is high |^ . Though for most systems careful applications of either procedure 
should yield results of comparable accuracy, they suffer from some common limitations. 
Firstly, elastic constants are obtained for a particular system size L. In order to take 
the infinite size (thermodynamic) limit one needs to simulate a sequence of systems with 
increasing values of L - often a computationally expensive proposition since equilibration of 
large systems take increasingly larger times. Secondly, these procedures are not general 
and fail for model systems of particles interacting via singular potentials eg., the hard 
sphere [0,0 or the hard disk [|18| system where the instantaneous force on a particle 



is not well defined. To obtain elastic constants of hard systems, one therefore needs to 



develop special methods [1(:,17]. For instance, the elastic constants of the hard sphere system 



was calculated by Runge and Chester |T^ by generalizing a technique used previously for 



calculating the hard sphere pressure |]T9|]. In this approach one tries to evaluate ensemble 
averages of quantities involving delta functions by calculating acceptance probabilities of 
virtual Monte Carlo steps. Such methods are cumbersome to use and may require ill defined 



averages of quantities whose variance has weak divergences [|T7|. In contrast, the method 



presented here has several advantages. Firstly, elastic constants are obtained from a coarse 
graining procedure which automatically, gives the infinite system values. Secondly, our 
procedure needs only the instantaneous particle configurations and makes no reference to 
the potential or forces. It is therefore quite generally applicable to any system for which 
these configurations or "snapshots" are known. Lastly, it will soon be evident that this 
technique is easy to use requiring a computational effort not much more than a calculation 
of, say, the pair distribution function for a given particle configuration. 
There are three essential elements or steps in the method, 

1. a procedure for calculating elastic strains from configurations. 

2. the coarse graining procedure of averaging fiuctuations of these strains over larger and 



larger sub -blocks of length L;, < L of a system of total size L and calculating the 
values of these fluctuations in the thermodynamic limit 

3. and, finally, converting the data for the strain fluctuations into elastic constants using 
well known results of continuum elasticity theory []7|,pO . 



We take up the second of these steps first as it is quite general and applicable to fluctua- 
tions of any intensive variable. Indeed, variants of this method have been used extensively 
2l|- p4t^ in the past for obtaining finite size scaled susceptibilities and goes by the name of 



"block - analysis". In the next section (section II) we describe this block -analysis technique 
for obtaining the susceptibility of the two -dimensional Ising model ||2^ for illustration. In 
section III we define the microscopic strains ejj(r) and describe how we obtain them from our 
computer simulations. We then explain how to adopt the coarse graining scheme described 
in section II for strains to finally obtain the elastic constants. This is followed by our results 
for the hard disk and the inverse twelveth power "soft" disk systems in two dimensions. We 
conclude this paper with a discussion of these results and enumerating future directions. 

II. FINITE SIZE SCALING OF FLUCTUATIONS IN SUBSYSTEMS 

The block analysis method has been used to obtain the compressibility of the Ising lattice 
gas |2^, the two -dimensional Lennard -Jones fluid [pSJ^ and fluids with internal classical 



and quantum degrees of freedom [p8| , p9| . Below, we describe a version of this method 
which allows us to explicitly and systematically incorporate finite size effects arising both 
from a non -zero correlation length as well as special constraints due to the nature of the 



ensemble used. This treatment is analogous to the one used in Ref. ^\ for analyzing density 
fluctuations in the two -dimensional Lennard Jones fluid. 

To begin, consider a general system described by a scalar order parameter 0(r). We are 
interested in obtaining thermodynamic properties of this system in the disordered phase. 
It is therefore sufficient to use the following quadratic Helmholtz free energy functional, 



F = kBT I d'^r Qr02 + lc{ V0(r)}2) (1) 

The correlation function G^^{q) implied by the above free energy in the high temperature 
phase (< >= 0) is given |§| by the following Ornstein -Zernicke form, 

f3G,4q) =< 0,0_, >= A'- ^ (2) 

Where ^^^ is the correlation length (= Jc/r), the (infinite system) susceptibility X^ =< 
0^ > is, in turn, given by \im.g^oPGtf,^{q), [3 = (ksT)^^ and the angular brackets < ... >= 
J[d<j>] exp(-/5F[{0}]). 

Consider further that we want to measure this susceptibility from a computer simulation 
of an Ising model within a finite simulation box of size L. In the course of the simulation 
we measure averaged within a sub-block of size Lb < L, 



L, 



'^ I ' d'^r(t){v). (3) 



Then the fluctuations of measured within this block is. 



< 02 >^^ L^f = U-" f" d^r'd'^r < 0(r)0(r') >, (4) 

= J"^' d'^rpC^^ir), 

Where G^(i,{r) is the inverse Fourier transform of the correlation function defined in Eq.(^ 
and is given p, quite generally, by, 

G,,{r)=r'X^^\r\('-'^Y{\r\/^^,) (5) 

where, 

pQQ p ^iz COS 

Y{r^)= z'^-'dz {2Ti)-''dn,——- (6) 

JO J [z^ + rj^\ 

One expects therefore that as, L^ ^ L the block susceptibility X.J^ -^ X?^. The be- 



haviour of the block susceptibilities, however, is strongly dependent on the ensemble [^ 
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which the simulation is carried out. For example, using the lattice gas pf language, in a 
"grand canonical" ensemble where the chemical composition cj) of the lattice gas is allowed 
to fluctuate keeping the chemical potential difference flxed (= in this case) the block sus- 
ceptibility does approach X^ for large L?,. In a canonical ensemble, however, the average of 
over the entire system is constrained to vanish for all times. In such a case the behaviour 
of X.i is more complicated but can, nevertheless, be explicitly determined as follows. 
Introducing the Lagrange multiplier k and defining the new correlator. 



G\r) = /'[#]</.(0)(/.(r)exp(-^F- /€ f d 



\ 



where the free energy F is given by Eq.(P), one shows that that G'^^{r) = G^ff,(r) — A^. 
The constant A^, is given by, 

AL = ^j'd'rG^^{r) (7) 

Repeating the analysis with G^(j,{r) replaced by G'^^{r) one obtains the desired finite size 
scaling form for X,^ in the canonical ensemble. We derive this form explicitly for the case 
oi d = 2 below. 

The correlation function G{r) is, 

G{r) = -r'X^Mr/O (8) 

TT 

where Kq is a Bessel function and we have suppressed the subscripts in G, X and ^ for 
simplicity. We have therefore, 

A^ = X^L'H{Lli) (9) 

with the function \I^(a) being defined as, 

^(a) = -a^ f [ dxdy KoiaJx^ + y^) (10) 

n Jo Jo * 

As can be easily verified \E'(a) goes to 1 within a range 0{1). Using the above expressions 
we find finally. 



p^L, ^ x'^[^>{xL/i) - ^{L/i)x\ (11) 

This is the main result of this section. 

Note that the form given above imphes X^'^ — > as Lf, ^ L. In general there will 



be higher order corrections to Eq.(ll) coming from higher order correlations signifying the 
breakdown of the quadratic form for the free energy Eq.(|l|). These introduce terms of order 
3;(2Af) implying extra parameters (with the constraint X^ = being imposed to eliminate 
one of them). In Fig. (1) we have plotted X^'' x x against x for various values of L/^. One 
observes that for L/C, ^ 1 ,^E'(a) -^ 1 and Eq.(|TT|) goes over to the following simple limiting 
form, 

X^^ X X = X°°[x - x^]. (12) 

One can then extract X°° from the slope of the linear region at small x. This is equivalent 
to the procedure followed for example in Ref. p4[ for extracting the compressibility of the 
lattice gas at the critical density. Obviously, for L/C, ~ 1, this construction becomes less 
well defined and the finding the "linear" region becomes subjective. Fitting the full data to 
the form given in Eq. ( pH]) makes it possible to extract X^ even when the system size is not 
much bigger than the correlation length. We illustrate this below using previously published 
data of Rovere, Nielaba and Binder |^ obtained from simulations of the lattice gas model 
using spin exchange dynamics (canonical ensemble). 

In Fig. (2) we plot X^'' as a function of L/Lb = 1/x. This choice of the axes is identical to 



the one used in the original work pj] and is equivalent. The points are the data of Ref. ||2J 
at various (dimensionless) temperatures T/J > T^/J = 2.2699. The curves through the data 
are fits to Eq.(|l2). The extracted compressibility X°° are plotted in Fig. (3) and compared 
with the previous estimates p3 based on linear extrapolation and the theoretical curve 



25| . Though, for large T where the correlation length C remains small one gets identical 
results, X°° extracted from the fits to the full curve continues to be closer to the theoretical 
(exact) estimate as T — ^ T^ (and consequently C, -^ oo). Simultaneously, we also obtain from 
Eq.(|TTl) an estimate for the correlation length C, which is shown in Fig. (4). The agreement 
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with theory ||26| is not as good. This is only to be expected since the correlation length is a 
more sensitive quantity than the susceptibility. Also, the fitting procedure cannot overcome 
errors built into the data because of critical slowing down near Tc which is known to lead to 



a systematic underestimation due to the use of a biased estimator pT |. 

Further, close to Tc the free energy functional (Eq. (|I])) ceases to be valid, the correlation 
function G(i)^{r) ~ \r\^i'^-'^+v) ^T^d this affects our estimate for the correlation length. This 
is particularly important in two dimensions |^ where, in the Ising lattice gas, the exponent 



r] = 1/4 indicates a rather large deviation of the correlation function from the Ornstein 
Zernike behavior, Eq. @), while in rf = 3, r^ ~ 0.03 and hence Eqs.(|^) and @ should be 
more accurate. 

Far away from a critical point, on the other hand, none of these criticisms apply and 
accurate estimates for the susceptibility can be obtained even from relatively small systems 
even m. d = 2 using the ideas described here. 

III. STRAIN FLUCTUATIONS AND ELASTIC CONSTANTS 

A. Theory 

Imagine a system in the constant NVT (canonical) ensemble at a fixed density p = 
N/V evolving in time t. For any "snapshot" of this system taken from this ensemble, the 
local instantaneous displacement field UR(t) defined over the set of lattice vectors {R} of a 
reference lattice (at the same density p) is 

Unit) = R(t) - R (13) 

where R(t) is the instantaneous position of the particle tagged by the reference lattice point 
R. In this paper we concentrate only on perfect crystalline lattices; if topological defects 
such as dislocations are present the analysis below needs to be modified. The instantaneous 
Lagrangian strain tensor etj defined at R is then given by P,^ , 



^'' ~ 2 \dR, dRi dRk dRj ) ^ ^ 

The strains considered here are always small and so we, hereafter, neglect the non -linear 
terms in the definition given above for simplicity. The derivatives are required at the ref- 
erence lattice points R and can be calculated by any suitable finite difference scheme once 
UR,(t) is known. We are now in a position to define coarse grained variables e^j' which are 
simply averages of the strain over a sub-block of size Lf,. The fluctuation of this variable 
then defines the size dependent compliance matrix Sijki =< eijtki >■ Before proceeding 
further, we introduce a compact Voigt notation (which replaces a pair of indices ij with one 
a) appropriate for two dimensional strains - the only case considered in this paper. Using 
1 = X and 2 = y we have, 

ij = 11 22 12 (15) 

a= 1 2 3 

The only nonzero components of the compliance matrix are 

Sll = < €xx^xx > = S22 (16) 

5*12 = < txx^yy > = S21 
'S'ss = 4 < exy^xy > 

It is also useful to define the following linear combinations 

5++ = <e+e+> =2{Sn + Su) (17) 

5__ = < e_e_ > = 2{Sn - ^12) 



Where e+ = exx + ^yy and e_ = exx — ^yy Once the block averaged strains e^j" are obtained, 
it is straight -forward to calculate these fiuctuations (for each value of Lf,). 

Since we are interested in the elastic properties of the system far away from any phase 
transition, a quadratic functional for the Helmholtz free energy F suffices. We therefore 
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use the following Landau functional appropriate for a two dimensional solid involving coarse 
grained strains to quadratic order in strains and its derivatives. 

F = Jd'^r {ciel + C2el + c^el + di{Ve+{r)f (18) 

+ci2(V6_(r))2 + d3(Ve3(r))n 

The coefficients q and di are, of course, not all independent because of further symmetries 
^J^ which are present for the two dimensional triangular solid though this does not influence 
the rest of our analysis. Note that, Eq. (|TBp is simply a sum of three independent functionals 
in e+, e_ and es each of the form given in Eq.([l|). Thus the analysis described in detail in 
section II carries over almost unchanged. There is one small difference, however. Even in 
the canonical ensemble with fixed box dimensions, the microscopic strain fiuctuations over 
the whole box are not zero but remains a small number of the order of (a/L)^ where a is 
the lattice parameter so that, 

j^^r < e,(r)e^(0) >= C^^p {^^\ (19) 

Incorporating this modification into Eq.(pJ]) we get. 

Si;;; = ^~ ^{xL/O - (^^L/O - C (I)') x' (20) 

+0(a;^). 

Where the index 7 takes the values +,— or 3 and x = L^/L and we have suppressed subscripts 
on C, and C for clarity. The above equation Eq.(pUD can now be used to obtain the system 
size independent quantities S*^, C, and C. 

Once the finite size scaled compliances are obtained the elastic constants viz. the Bulk 



modulus B = pdP/dp and the shear modulus /x are obtained simply using the formulae ||2^ 



(^^'=^-^P (22) 

/3^^=^-/3P (23) 
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where we assume that the system is under an uniform hydrostatic pressure P. The two 
expressions for fi should give identical results and constitutes an excellent internal check for 
numerical accuracy. 

This completes the description of our technique for obtaining elastic constants. It is 
obvious that throughout the derivation we do not ever refer to the interparticle interactions. 
The only place where this information is required is in the calculation of the pressure P of 
the system — a quantity routinely calculated in simulations. This, as we show later, is not 
a limitation even for systems with hard core potentials, where the calculation of pressure is 



more involved [jT^. We now describe our results for two model, two dimensional, solids. The 
elastic constants in each case are obtained for a high density perfect triangular solid. The 
generalization of this method for higher dimensions, more complicated (less symmetrical) 
crystal types and even for amorphous materials is straightforward. 

B. Hard disks 

Hard disks interact with the extremely short ranged potential, 

f3V{r) = for r > (T 
= cxD r < a, 

where a the hard sphere diameter can be used as the unit of length and there is no energy 
scale. This makes the hard disk system particularly easy to simulate and is a popular 
testing ground for theories and techniques 0. On the other hand, a calculation of elastic 
constants in this system is difficult since a harmonic description for such a solid does not 
exist at zero temperature. The energy vanishes and the free energy is wholly entropic in 
origin. For this reason most computational methods for obtaining elastic constants which 
work for smooth potentials have to be either discarded or modified in order to study elastic 
properties of hard systems [P"^-[T^ . The only previous study of the hard disk elastic moduli 
is by Wojciechowski and Brahkai []T3. These authors carried out a Monte Carlo simulation 
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of 56 hard disks in a box of variable shape in the constant stress ensemble [^ . The elastic 
constants were obtained from the fluctuations in the shape of the entire box. Results were 
checked for finite size effects by repeating a few test runs for 24, 30 and 90 particles. 

In our method the results derived in the previous sections carry over without any change 
to systems of particles interacting with hard potentials. Further we automatically calculate 
finite size scaled quantities. We present results for elastic constants of the hard disk system in 
two dimensions for a range of densities 0.95 < p* < 1.1 from Monte Carlo (MC) simulations 
with systems with sizes varying from 168 to 12480 particles. The simulations have been 
set up for a perfect triangular lattice in a slightly rectangular simulation box with periodic 
boundary conditions. The number of cells along each side of the box is adjusted to make the 
simulation box as close to a square as possible. In the hard disk system one can considerably 
accelerate the MC dynamics using special updating schemes 0. We use a square grid as an 
overlay on the simulation box, and choose the grid size to be small enough to accommodate 
only one disk. An occupancy list of the grid positions is computed and continuously updated. 
For an attempted move of a hard disk from one grid point to another, we first check if the 
new point is unoccupied and then check for overlap only with the particles occupying the 
neighboring grid points. At the density pa^ = p*= 1 a standard simulation had the length 
of 3x 10^ MC steps, every 10-th MC step has been taken for averaging observables, in 
particular the block analysis was done for 100 random placements of the block. 

In Fig. (5) we show the strain- strain fiuctuations as a function of the relative sub-block 
size computed in a system with A^ = 3120 hard disks at a density of p* = 1.0. The data are 
fitted to Eq.(^) where we keep terms up to order (L^/L)^. The fits to Eq.(pOD are excellent 
and the infinite system compliances S"^ can be determined immediately. To determine 
the sensitivity of our results to the total system size we have also simulated systems with 
A^ = 12480, 780 and an extremely small system consisting of A^ = 168 particles. Using the 
"infinite" system values determined from the A^ = 3120 particle system it is, in principle, 
possible to predict the behaviour of the strain fiuctuations S^^{L}j) for the other systems. In 
Fig. (6) we have plotted S^^^Lb) for 12480,3120,780 and 168 particles. The points are the 
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simulation data while the curves are fits to Eq. (pO[) where the value of 5*^ was constrained 
to be fixed at that obtained from the data for 3120 particles. The results are seen to be 
almost insensitive to the total system size as expected. 

Once the compliances 5*^^ are obtained the elastic constants, in units of ksT/a^, can be 
derived immediately using Eq. ( PT[ - P5D . Our results for densities other than p* = 1.0 were 
obtained from Monte Carlo simulations for A^ = 3120 hard disks. To obtain the pressure, 
which is required for evaluating the shear modulus /i, we simply integrate our bulk modulus 
B (independent of the pressure in two dimensions) starting from the rather high density 
p* = 1.1 where the free volume [^ expression for the pressure jSPa"^ = p* = 2p/(2/\^p—l) 



is accurate. Our results [Q for the equation of state for the hard disk system is shown in 
Fig. (7) and those for the elastic constants are shown in Fig. (8). The two expressions for 



the shear modulus in Eq. ( ]2^j23| ) give almost identical results and this gives us confidence 



about the internal consistency of our method. We have also compared our results to those 



of Wojciechowski and Brahka |T8[. We find that while their values of the pressure and 
bulk modulus are in good agreement with ours (and with free volume theory) they grossly 
overestimate the shear modulus. This is probably due to the extreme small size of their 
systems and/or insufficient averaging. Our results for the sub-block analysis shows that 
finite size effects are non -trivial for elastic strain fluctuations and they cannot be evaluated 
by varying the total size of the system from 24 to 90, an interval which is less than half of 
a decade |Q. One immediate consequence of our results is that the Cauchy relation p!8|j3l 



p = B/2 — P* is seen to be valid upto ±15% over the entire density range we studied (see 
Fig. (9) ) though there is a systematic deviation which changes sign going from negative 
for small densities to positive as the density is increased. This is in agreement with the 
usual situation in a variety of real systems [^ with central potentials and highly symmetric 
lattices and in disagreement with Ref. |]18|. We have also compared our estimates for the 
elastic constants with the density functional theory (DFT) of Rhysov and Tareyeva |^ . We 
flnd that both the bulk and the shear moduli are grossly overestimated - sometimes by as 
much as 100%. 
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C. Soft disks 

The system of particles in two dimensions interacting by a purely repulsive inverse 12*^^ 
power pair potential v{rij) of the form given by, 

v{r)=e(-j (24) 

has been studied [^^ quite extensively. This system has the advantage of being realistic 
without being too complicated, since the form of the potential ensures that the entire equa- 
tion of state can be determined from that of a single isotherm [^]. The quantities e and 
a sets the scales for energy and distance respectively and can be both set equal to unity. 
Both the zero and the finite temperature elastic constants of this system has been calculated 



over a large range of densities fS^]. We have used this system to test the applicablity of our 
method to molecular dynamics simulations. Our results here are not as extensive as in the 
hard disk case and we obtain elastic constants only for a single state point. 

We simulate this system with a simple leap -frog molecular dynamics code incorporating 
a Nose Hoover thermostat in order to generate configurations in the canonical NVT 
ensemble. The temperature T* = ksT/e is fixed at 1 and the density p* = pcr^ at 1.05 - a 
state sufficiently far from melting. The number of particles were chosen to be A^ = 780 which 
is same as that used in Ref. |3^ and corresponds to 26 x 30 unit cells of a triangular lattice 
within a nearly square box. Starting from the perfect lattice initial configuration, the system 
was equilibrated for more than 10^ molecular dynamics time steps At = 0.002 ( in units of 



ma'^/e where m is the mass of the particles). Subsequently, averages of thermodynamic 
quantities were calculated over ^ 10^ uncorrelated configurations. Our results for the elastic 
compliances S^^ are similar to that in the hard disk case and are shown in Fig. (10). 

The final estimates for the bulk modulus B = 77.96 and shear modulus /i = 23.34 (in 
units of kBT/a"^) compare well with those of Ref. |3^ (viz. B = 79.71 and /i = 24.96). 
Errors in our estimate for the bulk modulus arise from statistical error in the acquired data 
and from the fits this is around 3%. The shear modulus being a more sensitive quantity to 
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compute is less accurate and the two expressions for /i in Eq. (p2| , p3[) now differ by 10 — 15% 
(the number quoted above is the average value). This maybe due to the smaller size of our 
system compared to the hard disk case, and because subsequent configurations in molecular 
dynamics simulations are more correlated than in Monte Carlo. 

IV. CONCLUSION 

In summary, we show in this paper that a systematic coarse graining analysis of strain 
fluctuations yields elastic constants of solids from computer simulations to high accuracy. 
Our method incorporates finite size scaling and produces elastic constants in the thermo- 
dynamic limit. The procedure is simple to implement and is general enough to be carried 
out for any system without modification. Before we end this paper, the following comments 
are, perhaps, in order. 

The coarse graining variable:- Firstly, we introduced this work as a "test" case of a general 
coarse graining procedure which may be used to study the connection between microscopic 
computer simulations and long wavelength physics contained in continuum field theory ap- 
proaches for phase transitions in solids. In this regard we would like to point out that we 
found, as usual, the choice of the coarse graining variable to be important. Our results show 
that microscopic strains calculated by taking finite differences of the displacement field con- 
stitute the correct variable to coarse grain over. One could alternatively have averaged the 
displacement field u(r) and calculated the strains from the coarse grained u. This procedure 
happens to produce wrong results giving us elastic constants which are orders of magnitude 
too large. 

The strain correlation length:- One of the new results from our calculations is the correlation 
length ^a/B of strain fluctuations. This is found to be small - 2 — 3 lattice spacings - for 
all components of the strain-strain correlation functions and for both the systems. We 
have checked this independently by explicitly measuring the correlation function GsadlRH) 
defined for all the lattice vectors R of the two dimensional triangular reference lattice in the 
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soft sphere system. Though, the simple Ornstein Zernicke form is inadequate to describe 
detailed features of the correlation function and actual values for the correlation length is 
hard to estimate, preliminary results from our simulations do support the above contention. 
The correlation of the local density p(r) or its phase - the displacement field u(R) is, of 
course, long ranged, decaying algebraically as it should in the solid state. 
Renormalization by defects:- Our results for the elastic constants are obtained for high 
density perfect solids. In general a solid contains point (vacancies and interstitials) and 
line (dislocations) defects. For example, in the hard disk case dislocations start appearing 
in our systems below a density of p* = 0.95 a range we have not explored in this paper. 
In principle, there is no reason why our method cannot be adapted for systems containing 
defects though this involves considerable computational complexity. There are basically two 
problems which arise when one wishes to calculate elastic constants in the presence of defects. 
Firstly, we have to ensure that the density of each type of defect on the average attain their 
equilibrium value. This is nontrivial because nucleation barriers for defects densities are 
usually high which means large system sizes and long simulation times are required. Defect 
mobilities are sluggish in a solid but to ensure that they are fully equilibrated one has to 
wait long enough to allow a typical dislocation to travel a distance equal to the system size 
Il39| . Secondly, once we are sure that our configurations contain, on the average, the required 
defect densities, we have to evaluate the strain field in presence of these defects. This is, 
by far, the easier part. Defects can be either isolated point defects which arise in grand 
canonical simulations or vacancy -interstitial pairs which nucleate pairs of dislocations of 
opposite sign. The concentration of point defects in solids is vanishingly small ~ 10"^ — 10~^ 
atomic percent and they are not expected to change the elastic constants substantially. In 
any case they can be easily taken into account by redefining the reference lattice points. 
Topological defects introduce, in addition, a singular part in the displacement field so that 
the strain field cannot now be evaluated simply by taking numerical derivatives. However, 
the singular contribution for each dislocation is known analytically - so that for a given 
configuration one can locate topological defects of each type and treat the smooth and the 
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singular parts separately. Lastly, we would like to point out that the above two problems viz. 
(1) obtaining an equilibrium defect concentration and (2) evaluating the finite size scaling 
of singular strain fields of configurations containing defects are present in all techniques for 
calculating elastic constants although they seem not to have recieved so far the attention 
they deserve. 

The renormalization of elastic constants by dislocations can also be obtained approxi- 
mately by using standard recursion relations [jl2| once the core energy of a dislocation is 
known (for instance from a separate simulation) and the present (manifestly "bare" ) values 
of the elastic constants are used as inputs. Such a calculation has applications in a study 



of two dimensional melting |TT|JT^ . We are, at present, carrying out detailed calculations 
of the elastic constants, equation of state and dislocation core energies of the hard disk and 
inverse power triangular solids to investigate their melting behaviour. Further, near to the 
melting transition additional finite size effects would be manifest (due to diverging corre- 
lation lengths) and they have to be taken into account separately which would necessarily 
involve simulations of a much larger scale than employed here. For example, it is known 
that near the liquid to hexatic transition (p* = 0.89) A^ = 256^ hard disks are required 



in order to reach the thermodynamic limit. 

Evaluation of local stresses:- Lastly, we would like to point out that our procedure can 
also be used in an "inverse" mode where knowing the elastic constants for a system, strain 
fluctuations can be used to calculate local stresses. This is especially useful in experimental 



studies |41,42] of melting behaviour of colloidal particles where high quality digitized particle 



images are available. Efforts in this direction are in progress. 
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FIGURES 
FIG. 1. A plot of X^'' as a function of the relative sub-block size Li,/L as given by Eq.(ll) 
for various values of the parameter L/^ (= 1, 5, 10 &: 50). 



FIG. 2. The susceptibility A!^''{T) of a two dimensional Ising system in the constant magne- 
tization ensemble as a function of the relative sub-block size L/Li, for values of the temperature 
T/J > Tc/J. The symbols refer to simulation data of the critical lattice gas simulation of Ref. (23) 
and the lines are fits of Eq. ([Tl|) to these data. 



FIG. 3. The estimates for XooiT) as obtained from the fits shown in Fig. 2. (O) compared 
with the estimates of Ref. (23) (+). Note that the finite size scaling analysis described here yields 
estimates which are closer to the exact result ( curve ). 



FIG. 4. The correlation length ^ (in units of the lattice parameter a) as a function of the 
scaled temperature T/J of the lattice gas model at the critical density. Symbols: our estimates 



from fits to the data of Ref. |24] and curve: theory Ref. pq]. 



FIG. 5. Strain strain fluctuations or elements of the compliance matrix S^yy defined in Eq. 
([T6|) (7 = +, — or 3) as a function of relative sub-block size Lb/L, symbols: simulation data; curve: 
fit to the scaling function Eq. (pC|). The results (symbols) shown are for a system of iV = 3120 
hard disks at p* = 1.0. 
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FIG. 6. The infinite system susceptibility S'33 obtained from the data for N = 3120 particles 
is used to predict the finite size behaviour of A^ = 168, 780 and 12480 particles. For A^ = 168, 780 
and 3120 the symbols are simulation data and the solid lines are the best fits to the form given in 
Eq. (|20| ) where the parameter S'33 i^ kept fixed at the value obtained from fits to the data for 3120 
particles. For A^ = 12480 we have acquired data (symbols) only for L/Lj = an integer (between 4 
and 35) and the dotted line is a straight line with the slope given by the same S'33. 

FIG. 7. The equation of state of the hard disk solid, pressure P* as a function of the density 
p*. We compare our results (O) with those of Ref. ||l^ (+) and free volume theory (line). 



FIG. 8. The bulk {B) and shear (fi) moduli in units of ksT/a'^ for the hard disk solid. Our 
results for B (/i) are given by □ (O). The values for the corresponding quantities from Ref. [^] are 
given by + and x . The line through the bulk modulus values is the analytical expression obtained 
from the free volume prediction for the pressure. The line through our shear modulus values is 
obtained from the free volume bulk modulus using the Cauchy relation ^u = B/2 — P. The error 
bars in our values for the shear modulus correspond to the two alternative formulae for evaluating 



fi as given in Eq.(22,23) and represents our most conservative error estimates. 



FIG. 9. The percentage deviation Ac of our shear modulus values /x from that obtained from 
our bulk moduli B using the Cauchy relation fi = B/2 — Pas a funtion of the density p*. The error 
bars in this graph correspond to the two formulae for evaluating p as given in Eq.( P2| , P3 ) a-s in the 
previous figure. 

FIG. 10. Same as in Fig. 5 but for N = 780 soft disks interacting with the inverse 12*'* power 
potential at p* = 1.05 and T* = 1. 
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